%%%% Solution Algorithm for the Representative-Agent Equivalent for the Aiyagari (1994,QJE) Model for Welfare Calculations 
%%%% Written by Orhan Torul & Oguz Oztunali
%%%% October 2017, Bogazici University, Istanbul

%%
%%
%% 0. housekeeping
clc;                                % clear screen
clear;                              % clear memory
close all;                          % close open windows
%%
%% 1. parameters and functional forms

% parameters
alpha  = 0.56;                              % capital income share from Penn World Table 9.0 (PWT)
b      = 0;                                 % borrowing constraint as in Aiyagari (1994, QJE)
beta   = 89/100;                            % subjective discount factor from PWT
delta  = 5.5/100;                           % depreciation rate of physical capital from PWT
gamma  = 3/2;                               % inverse elasticity of intertemporal substitution from macroeconomics literature
varphi = 2/3;                               % Frisch elasticity of labor supply from macroeconomics literature
zbar   = 0.783472854785600;                 % effective labor productivity from the heterogeneous-agent incomplete market model
utilityaiyagari= -1.996968282436092;        % average utilitarian contemporaneous utility in the heterogeneous-agent environment
utilityaiyagariss= utilityaiyagari/(1-beta);% infinite summation of utilityaiyagari with subjective discounting     


% utility and (inverse) marginal utility functions
% contemporaneos utility is defined as u(c)-v(x) with the functinal forms below

u     = @(c) (c.^(1-gamma))/(1-gamma);          % utility from consumption
v     = @(h) (h.^(1+1/varphi))/(1+1/varphi);    % disutility from labor supply
up    = @(c) c.^(-gamma);                       % marginal utility of consumption
invup = @(x) x.^(-1/gamma);                     % inverse of marginal utility of consumption
vp    = @(h) h.^(1/varphi);                     % marginal disutility of labor supply
invvp = @(x) x.^(varphi);                       % inverse marginal disutility of labor supply 



hbar = 0.912;                  % initial guess for hours worked, which delivers the same average utilitarian contemporaneous utility as in the heterogeneous-agent environment
utilitydifference= 10;
crit             = 0.0001;

while utilitydifference>crit
%% 2. deterministic steady-state
 
% initial guess
y0 = [10, 0.3]; 
% assignment of endogenous variables are as follows:
% y(1) capital (k)
% y(2) forgone lump-sum consumption (closs)

% option setting for the fsolve function
options= optimset('maxiter',999999999999','tolfun',10e-8,'tolx',10e-8);
[SS, fval,exitflag,output]=fsolve(@(y) steady_counterfactual(y,alpha,beta,delta,gamma,varphi,zbar,hbar),y0,options);

kss=SS(1);
closs=SS(2);
hss=hbar;
lss=hss*zbar;

css=((hss^(1/varphi))/((1-alpha)*hss^(-alpha)*(kss^alpha)*zbar^(1-alpha)))^(-1/gamma);
wss=(1-alpha)*kss^alpha*lss^(-alpha);
rss=alpha*kss^(alpha-1)*lss^(1-alpha)-delta;
gdpss=kss^alpha*lss^(1-alpha);
incomess=gdpss-delta*kss;

disp ('steady-state for the decentralized economy:')
disp('-------------------------------------------------')
fprintf('output (y):              %.4f \n',gdpss);
fprintf('capital (k):             %.4f \n',kss);
fprintf('labor (h):               %.4f \n',hss);
fprintf('effective labor (l):     %.4f \n',lss);
fprintf('consumption (c):         %.4f \n',css);
fprintf('wage (w):                %.4f \n',wss);
fprintf('real interest rate (r):  %.4f \n',rss);
disp('-------------------------------------------------')
fprintf('utility from c:          %.4f \n',u(css));
fprintf('disutility from n:       %.4f \n',v(hss));
fprintf('contemporaneous utility: %.10f \n',u(css)-v(hss));

utilitycounterfactual=u(css)-v(hss);
utilitycounterfactualss=utilitycounterfactual/(1-beta);
utilitydifference=abs(utilitycounterfactualss-utilityaiyagariss)

hbar=hbar+0.00001  %updating guess in the upward direction
end

% exporting to LaTeX
meanmatrix=[kss hss zbar lss gdpss css rss wss rss*kss wss*lss ];
Det_TableStationaryMeans=[meanmatrix];
columnlabelStationaryMeans = {'$K$', '$H$' ,'$Z$' ,'$L$', '$Y$', '$C$', '$r$', '$w$' ,'$rK$', '$wL$'};
matrix2latex(Det_TableStationaryMeans, 'Det_TableStationaryMeansCounterfactual.tex', 'columnLabels', columnlabelStationaryMeans,'format', '%-6.3f', 'size', 'footnotesize');


closspercent=closs/(css+closs);              % forgone consumption (as a percentage of steady-state consumption) for indifference

cssbenchmark = 1.897274725917953;            % steady-state consumption in the benchmark representative-agent model with endogenous labor supply

hssbenchmark= 0.683342992989747;             % steady-state hours worked in the benchmark representative-agent model with endogenous labor supply

c_relative=1-css/cssbenchmark;               % forgone steady-state consumption relative to the benchmark model environment
h_relative=hss/hssbenchmark-1;               % extra hours worked relative to the benchmark model environment           
